topo_change <- st_read('data/topo_change/topo_change_polygons.shp')
## Reading layer `topo_change_polygons' from data source `C:\Users\jesca\OneDrive\Documents\WR674\mining_impacts_to_elevation\data\topo_change\topo_change_polygons.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8958 features and 51 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -2332814 ymin: 407783 xmax: 2238422 ymax: 3128340
## epsg (SRID): NA
## proj4string: +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
#Checkout the column names
#names(topo_change)
biggest_area <- topo_change %>%
arrange(desc(AREA_SQ_KM)) %>%
slice(1:10)
# Map all mines
mapview(biggest_area)
# Subset to just Arizona
az_mine <- biggest_area %>%
filter(QUADNAME == 'esperanza_mill_AZ')
#Check that it is the right site
mapview(az_mine)
#Check projection of az_mine
#st_crs(az_mine)
az_raster_before <- get_elev_raster(az_mine,z=12)
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
#Look at the structure of the data
#str(az_raster)
#Summary of the data
#summary(az_raster)
plot(az_raster_before)
#mapview(az_raster_before)
small_az <- aggregate(az_raster_before,fact=3)
mapview(small_az)
az_slope <- terrain(small_az,opt='slope',unit='degrees')
mapview(az_slope)
az_aspect <- terrain(small_az,opt='aspect',unit='degrees')
mapview(az_aspect)
az_matrix <- matrix(raster::extract(small_az,
raster::extent(small_az),buffer=500),
nrow=ncol(small_az), ncol=nrow(small_az))
az_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(az_matrix),color='desert') %>%
add_shadow(ray_shade(az_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(az_matrix)) %>%
plot_3d(az_matrix,zscale=10,fov = 0,theta=135,zoom=0.75)
render_snapshot(clear = TRUE)
For this section we are really going to focus on the mine area so first I’m going to buffer our az_mine polygon by 2km and then clip both the before and after mining rasters to this new shape
#Read in raster data
az_raster_after <- raster('data/srtm_14_06.tif')
az_mine_2km <- az_mine %>%
st_transform(2163) %>%
st_buffer(2000) %>%
#Reproject into SRTM datum (WGS84)
st_transform(st_crs(az_raster_after))
#Get extent of buffered mine area
aoi <- extent(az_mine_2km)
#Crop after to area of interest
mine_after <- crop(az_raster_after,aoi)
#Crop before to area of interest (need to also reproject)
mine_before <- crop(az_raster_before %>%
projectRaster(.,crs=projection(mine_after)),
aoi) %>%
projectRaster(.,mine_after) #match resolution
#Force rasters to be exactly same size
mine_before_same <- trim(mine_before)
mine_after_same <- crop(mine_after,mine_before_same)
#Stack these things in a rasterBrick (which we can do because we
#made them the exact same size)
mines <- brick(mine_before_same,mine_after_same)
#Add a dem_difference layer
mines[[3]] <- mines[[1]] - mines[[2]]
#Rename raster layers
names(mines) <- c('Pre_mining','Post_mining','DoD')
plot(mines)
All data for this section should be the data clipped to the mining area (generated in the code chunk above)
after_slope <- terrain(mine_after_same,opt='slope',unit='degrees')
plot(after_slope)
before_slope <- terrain(mines[[1]],opt = 'slope', unit = 'degrees')
plot(before_slope)
diff_slope <- terrain(mines[[3]], opt = 'slope', unit = 'degrees')
#plot(diff_slope)
mine_slopes <-brick(before_slope, after_slope, diff_slope)
names(mine_slopes) <- c('Pre_mining','Post_mining','DoD')
plot(mine_slopes)
pre_matrix <- matrix(raster::extract(mine_before_same,
raster::extent(mine_before_same),buffer=500),
nrow=ncol(mine_before_same), ncol=nrow(mine_before_same))
pre_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(pre_matrix),color='desert') %>%
add_shadow(ray_shade(pre_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(pre_matrix)) %>%
plot_3d(pre_matrix,zscale=10,fov = 0,theta=135,zoom=0.85)
render_snapshot(clear = TRUE)
post_matrix <- matrix(raster::extract(mine_after_same,
raster::extent(mine_after_same),buffer=500),
nrow=ncol(mine_after_same), ncol=nrow(mine_after_same))
dim(post_matrix)
## [1] 88 80
dim(pre_matrix)
## [1] 88 80
post_matrix %>%
sphere_shade(texture = 'desert') %>%
add_water(detect_water(post_matrix),color='desert') %>%
add_shadow(ray_shade(post_matrix,zscale=2),0.5) %>%
add_shadow(ambient_shade(post_matrix)) %>%
plot_3d(post_matrix,zscale=10,fov = 0,theta=135,zoom=0.85)
render_snapshot(clear = TRUE)